This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
library(raster)
size_dat=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/fire_growth_5days_v4.txt", header=T, row.names=NULL)
size_dat=as.data.frame(size_dat)
size_dat=size_dat[-1,]
size_dat=size_dat[,-1]
colnames(size_dat)=c("firename","year","cause","size1","size2","size3","size4","size5","final_firesize","mean_precip1","mean_precip2","mean_precip3","mean_precip4","mean_precip5","mean_tmax1","mean_tmax2","mean_tmax3","mean_tmax4","mean_tmax5","mean_tmean1","mean_tmean2","mean_tmean3","mean_tmean4","mean_tmean5","mean_vpdmax1","mean_vpdmax2","mean_vpdmax3","mean_vpdmax4","mean_vpdmax5","mean_windspeed1","mean_windspeed2","mean_windspeed3","mean_windspeed4","mean_windspeed5","landcover","ecosystem","biomass","elevation")
size_dat2 <- data.frame(lapply(size_dat, function(x) as.numeric(as.character(x))))
NAs introduced by coercion
size_dat2$human = 0
size_dat2$human[size_dat2$cause !=1 & size_dat2$cause !=14 & size_dat2$cause !=17]=1
size_dat2$human[size_dat2$cause ==1 ]=2
pro1 =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1 ), ]
pro2 =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1 ), ]
t.test(pro1$size1,pro2$size1)
Welch Two Sample t-test
data: pro1$size1 and pro2$size1
t = -2.0397, df = 74.381, p-value = 0.04493
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-22.0986230 -0.2593373
sample estimates:
mean of x mean of y
3.056385 14.235365
t.test(pro1$size2,pro2$size2)
Welch Two Sample t-test
data: pro1$size2 and pro2$size2
t = -2.758, df = 73.138, p-value = 0.007341
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-32.702045 -5.266208
sample estimates:
mean of x mean of y
8.254394 27.238521
t.test(pro1$size3,pro2$size3)
Welch Two Sample t-test
data: pro1$size3 and pro2$size3
t = -2.9002, df = 58.203, p-value = 0.005254
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-51.965993 -9.526741
sample estimates:
mean of x mean of y
14.10484 44.85121
t.test(pro1$size4,pro2$size4)
Welch Two Sample t-test
data: pro1$size4 and pro2$size4
t = -3.1797, df = 53.078, p-value = 0.002461
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-74.99798 -16.98040
sample estimates:
mean of x mean of y
17.68430 63.67349
t.test(pro1$size5,pro2$size5)
Welch Two Sample t-test
data: pro1$size5 and pro2$size5
t = -3.4218, df = 39.443, p-value = 0.001462
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-117.71273 -30.26844
sample estimates:
mean of x mean of y
19.71795 93.70853
pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1), ]
length(pro$year)
[1] 68
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1), ]
length(pro$year)
[1] 77
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 2), ]
length(pro$year)
[1] 18
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 2), ]
length(pro$year)
[1] 9
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
length(!is.na(pro$size5))
[1] 9
pro =size_dat2[which(size_dat2$human == 1 & size_dat2$ecosystem == 6), ]
length(pro$year)
[1] 41
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 2& size_dat2$ecosystem == 6), ]
length(pro$year)
[1] 79
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
NA
NA
plot fire size map QGIS
library(data.table)
DT= data.table(res)
fire_size = DT[ , .SD[which.min(growth_km)], by = firename]
fire_size = fire_size[,c()]
spread_list =list.files(viirs_dir, pattern = "_daily.shp$", recursive = TRUE, full.names=T)
i=1
l=shapefile(spread_list[i])
for (i in 2:length(spread_list)){
p=shapefile(spread_list[i])
l <- rbind(l, p, makeUniqueIDs = TRUE)
}
l$human = 0
l$human[l$cause !=1 & l$cause !=14 & l$cause !=17]=1
l$human[l$cause ==1 ]=2
writeOGR(l, "/Users/stijnhantson/Documents/projects/VIIRS_ros/", layer= "all_fires", driver="ESRI Shapefile", overwrite_layer = T)
library(raster)
dr =shapefile("/Users/stijnhantson/Documents/projects/VIIRS_ros/frap_subset.shp")
dr1 =shapefile("/Users/stijnhantson/Documents/data/FRAP/fire18_1.shp")
dr1$YEAR_=as.numeric(as.character(dr1$YEAR_))
dr1$Shape_Area=as.numeric(as.character(dr1$Shape_Area))
dr1=dr1[!is.na(dr1$YEAR_), ]
dr1=dr1[dr1$YEAR_>2011,]
sum(dr1$Shape_Area)
[1] 25651974971
sum(dr$Shape_Area)
[1] 22801510990
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/mean_vs_max_ros_v1.tif", width = 5, height = 5, units = 'in', res = 300)
plot(res$ros_km,res$ros_mean_km, xlim=c(0,25),ylim=c(0,10), xlab="maximum fire-spread-rate (km/day)",ylab="mean fire-spread-rate (km/day)", cex.lab=1.4,cex.axis = 1.3)
dev.off()
null device
1
plot mean vs max fire rate-of-spread
#summary(res)
plot(res$ros_km,res$ros_mean_km, xlab="maximum fire-spread-rate (km/day",ylab="mean fire-spread-rate (km/day)")
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/mean_vs_max_ros_v1.tif", width = 5, height = 5, units = 'in', res = 300)
plot(res$ros_km,res$ros_mean_km, xlim=c(0,25),ylim=c(0,10), xlab="maximum fire rate-of-spread (km/day)",ylab="mean fire rate-of-spread (km/day)", cex.lab=1.3,cex.axis = 1.25)
dev.off()
quartz_off_screen
2
difference between human and lightnign fires
me=0
me1=0
days = c("day1","day2","day3","day4","day5")
pro1 = res[res$fire_day == 1 & res$human == 1,]
pro2 = res[res$fire_day == 2 & res$human == 1,]
pro3 = res[res$fire_day == 3 & res$human == 1,]
pro4 = res[res$fire_day == 4 & res$human == 1,]
pro5 = res[res$fire_day == 5 & res$human == 1,]
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
pro1h = res[res$fire_day == 1 & res$human == 2,]
pro2h = res[res$fire_day == 2 & res$human == 2,]
pro3h = res[res$fire_day == 3 & res$human == 2,]
pro4h = res[res$fire_day == 4 & res$human == 2,]
pro5h = res[res$fire_day == 5 & res$human == 2,]
me1[1] =mean(pro1h$growth_km,na.omit=T)
me1[2] =mean(pro2h$growth_km,na.omit=T)
me1[3] =mean(pro3h$growth_km,na.omit=T)
me1[4] =mean(pro4h$growth_km,na.omit=T)
me1[5] =mean(pro5h$growth_km,na.omit=T)
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure1_v1.tif", width = 10, height = 5, units = 'in', res = 300)
par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^-1*')'),ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)
dev.off()
quartz_off_screen
2
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_figure1_v1.tif", width = 10, height = 5, units = 'in', res = 300)
par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^-1*')'),ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)
dev.off()
quartz_off_screen
2
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
t.test(log(pro1$growth_km),log(pro1h$growth_km))
Welch Two Sample t-test
data: log(pro1$growth_km) and log(pro1h$growth_km)
t = 6.5083, df = 169.36, p-value = 8.265e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
1.492386 2.791866
sample estimates:
mean of x mean of y
1.3226886 -0.8194377
t.test(log(pro2$growth_km),log(pro2h$growth_km))
Welch Two Sample t-test
data: log(pro2$growth_km) and log(pro2h$growth_km)
t = 6.9843, df = 134.4, p-value = 1.199e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
1.369607 2.451706
sample estimates:
mean of x mean of y
2.7559564 0.8453003
t.test(log(pro3$growth_km),log(pro3h$growth_km))
Welch Two Sample t-test
data: log(pro3$growth_km) and log(pro3h$growth_km)
t = 5.3694, df = 137.17, p-value = 3.286e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.9293526 2.0129234
sample estimates:
mean of x mean of y
3.186372 1.715234
t.test(log(pro4$growth_km),log(pro4h$growth_km))
Welch Two Sample t-test
data: log(pro4$growth_km) and log(pro4h$growth_km)
t = 4.6217, df = 122.92, p-value = 9.467e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.6909157 1.7261008
sample estimates:
mean of x mean of y
3.513204 2.304696
t.test(log(pro5$growth_km),log(pro5h$growth_km))
Welch Two Sample t-test
data: log(pro5$growth_km) and log(pro5h$growth_km)
t = 4.9328, df = 103.98, p-value = 3.089e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.9170094 2.1499736
sample estimates:
mean of x mean of y
3.895266 2.361774
##################3 for western cordillera ecoregion ##################
me=0
me1=0
pro1 = res[res$fire_day == 1 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)
plot(me,type="o", ylim=c(0,150),ylab="fire size (km2)",xlab="", xaxt='n')
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o")
t.test(log(pro1$growth_km),log(pro1h$growth_km))
Welch Two Sample t-test
data: log(pro1$growth_km) and log(pro1h$growth_km)
t = 3.2086, df = 69.426, p-value = 0.002019
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.5363985 2.2992588
sample estimates:
mean of x mean of y
0.4219047 -0.9959240
t.test(log(pro2$growth_km),log(pro2h$growth_km))
Welch Two Sample t-test
data: log(pro2$growth_km) and log(pro2h$growth_km)
t = 4.9267, df = 97.316, p-value = 3.427e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.9938201 2.3346153
sample estimates:
mean of x mean of y
2.303104 0.638886
t.test(log(pro3$growth_km),log(pro3h$growth_km))
Welch Two Sample t-test
data: log(pro3$growth_km) and log(pro3h$growth_km)
t = 3.4468, df = 82.843, p-value = 0.0008936
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.4842092 1.8055445
sample estimates:
mean of x mean of y
2.695279 1.550402
t.test(log(pro4$growth_km),log(pro4h$growth_km))
Welch Two Sample t-test
data: log(pro4$growth_km) and log(pro4h$growth_km)
t = 3.019, df = 73.938, p-value = 0.003479
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.3104149 1.5156656
sample estimates:
mean of x mean of y
3.092997 2.179957
t.test(log(pro5$growth_km),log(pro5h$growth_km))
Welch Two Sample t-test
data: log(pro5$growth_km) and log(pro5h$growth_km)
t = 3.1504, df = 55.57, p-value = 0.002625
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.4138127 1.8597722
sample estimates:
mean of x mean of y
3.402659 2.265867
mean(pro1$growth_km,na.omit=T)
[1] 16.60532
mean(pro1h$growth_km,na.rm=T)
[1] 2.714784
for mediteranean california
pro1 = res[res$fire_day == 1 & res$human == 1 & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & (res$eco1 == 11),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & (res$eco1 == 11),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
plot(me,type="o", ylim=c(0,150),ylab="fire size (km2)",xlab="", xaxt='n')
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o")
for difference in autumn
pro1 = res[res$fire_day == 1 & res$human == 1 & res$month >8,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month >8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month >8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month >8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month >8,]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & res$month >8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month >8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month >8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month >8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month >8,]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
plot(me,type="o", ylim=c(0,150),ylab="fire size (km2)",xlab="", xaxt='n')
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o")
for difference in summer
pro1 = res[res$fire_day == 1 & res$human == 1 & res$month <=8 ,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month <=8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month <=8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month <=8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month <=8,]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & res$month <=8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month <=8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month <=8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month <=8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month <=8,]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
plot(me,type="o", ylim=c(0,150),ylab="fire size (km2)",xlab="", xaxt='n')
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o")
NA
NA
for difference in summer in western cordillera
pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
t.test(log(pro1$growth_km),log(pro1h$growth_km))
Welch Two Sample t-test
data: log(pro1$growth_km) and log(pro1h$growth_km)
t = 1.3062, df = 51.551, p-value = 0.1973
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.3236595 1.5299715
sample estimates:
mean of x mean of y
-0.1773847 -0.7805406
t.test(log(pro2$growth_km),log(pro2h$growth_km))
Welch Two Sample t-test
data: log(pro2$growth_km) and log(pro2h$growth_km)
t = 3.6535, df = 78.926, p-value = 0.0004639
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.5639969 1.9140943
sample estimates:
mean of x mean of y
2.0787175 0.8396718
t.test(log(pro3$growth_km),log(pro3h$growth_km))
Welch Two Sample t-test
data: log(pro3$growth_km) and log(pro3h$growth_km)
t = 2.1785, df = 57.055, p-value = 0.03352
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.05727942 1.36024900
sample estimates:
mean of x mean of y
2.551723 1.842958
t.test(log(pro4$growth_km),log(pro4h$growth_km))
Welch Two Sample t-test
data: log(pro4$growth_km) and log(pro4h$growth_km)
t = 2.353, df = 57.21, p-value = 0.02208
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1074894 1.3347908
sample estimates:
mean of x mean of y
2.992332 2.271192
t.test(log(pro5$growth_km),log(pro5h$growth_km))
Welch Two Sample t-test
data: log(pro5$growth_km) and log(pro5h$growth_km)
t = 2.3446, df = 41.508, p-value = 0.02391
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.105315 1.410281
sample estimates:
mean of x mean of y
3.313679 2.555881
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
plot(me,type="o", ylim=c(0,150),ylab="fire size (km2)",xlab="", xaxt='n')
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o")
pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
`
how many days does it take to reach 75% burnt area
res75 = res[res$per_ba > 0.75,]
peak_day = as.data.frame(aggregate(res75$fire_day, by = list(res75$firename,res75$cause), min))
peak_day=subset(peak_day,x < 55)
hi = hist(peak_day$x,prob =F, breaks= c(0:53), xlim=c(0,55), ylab="number of fires", xlab="days", cex.lab=1.4,cex.axis=1.3)
out1 = subset(res75,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 )
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
quantile(peak_day1$x,0.50,type=3)
50%
10
quantile(peak_day2$x,0.50,type=3)
50%
3
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
tiff(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/time_to_reach75_v1.tif",width=2500,height=2000, res=350)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
3
difference in fire size for first 5 days across california and both ecosystems
res$ros1 = res$max_ros+1
out1 = subset(res,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )
hum1 = out2[out2$fire_day ==1,]
hum2 = out2[out2$fire_day ==2,]
hum3 = out2[out2$fire_day ==3,]
hum4 = out2[out2$fire_day ==4,]
hum5 = out2[out2$fire_day ==5,]
lig1 = out1[out1$fire_day ==1,]
lig2 = out1[out1$fire_day ==2,]
lig3 = out1[out1$fire_day ==3,]
lig4 = out1[out1$fire_day ==4,]
lig5 = out1[out1$fire_day ==5,]
mean(hum1$growth)
[1] 18753898
mean(hum2$growth)
[1] 36707455
mean(hum3$growth)
[1] 57176147
mean(hum4$growth)
[1] 80694648
mean(hum5$growth)
[1] 115079142
mean(lig1$growth)
[1] 2875395
mean(lig2$growth)
[1] 10065897
mean(lig3$growth)
[1] 18876394
mean(lig4$growth)
[1] 28689523
mean(lig5$growth)
[1] 37533565
t.test(log10(hum1$growth),log10(lig1$growth))
Welch Two Sample t-test
data: log10(hum1$growth) and log10(lig1$growth)
t = 6.5083, df = 169.36, p-value = 8.265e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.6481351 1.2124922
sample estimates:
mean of x mean of y
6.574436 5.644123
t.test(log10(hum2$growth),log10(lig2$growth))
Welch Two Sample t-test
data: log10(hum2$growth) and log10(lig2$growth)
t = 6.9843, df = 134.4, p-value = 1.199e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.5948126 1.0647622
sample estimates:
mean of x mean of y
7.196897 6.367109
t.test(log10(hum3$growth),log10(lig3$growth))
Welch Two Sample t-test
data: log10(hum3$growth) and log10(lig3$growth)
t = 5.3694, df = 137.17, p-value = 3.286e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.4036127 0.8742015
sample estimates:
mean of x mean of y
7.383824 6.744917
t.test(log10(hum4$growth),log10(lig4$growth))
Welch Two Sample t-test
data: log10(hum4$growth) and log10(lig4$growth)
t = 4.6217, df = 122.92, p-value = 9.467e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.3000609 0.7496361
sample estimates:
mean of x mean of y
7.525765 7.000917
t.test(log10(hum5$growth),log10(lig5$growth))
Welch Two Sample t-test
data: log10(hum5$growth) and log10(lig5$growth)
t = 4.9328, df = 103.98, p-value = 3.089e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.3982521 0.9337217
sample estimates:
mean of x mean of y
7.691692 7.025706
hum1 = out2[out2$fire_day ==1 & out2$eco1==6,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==6,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==6,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==6,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==6,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==6,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==6,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==6,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==6,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==6,]
mean(hum1$growth)
[1] 16605316
mean(hum2$growth)
[1] 30366275
mean(hum3$growth)
[1] 44666667
mean(hum4$growth)
[1] 63374279
mean(hum5$growth)
[1] 85629090
mean(lig1$growth, na.rm=T)
[1] 2747161
mean(lig2$growth, na.rm=T)
[1] 8686182
mean(lig3$growth, na.rm=T)
[1] 15363041
mean(lig4$growth, na.rm=T)
[1] 21997987
mean(lig5$growth, na.rm=T)
[1] 26941885
t.test(log10(hum1$growth),log10(lig1$growth))
Welch Two Sample t-test
data: log10(hum1$growth) and log10(lig1$growth)
t = 3.1605, df = 69.967, p-value = 0.002328
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.2244363 0.9921939
sample estimates:
mean of x mean of y
6.183231 5.574916
t.test(log10(hum2$growth),log10(lig2$growth))
Welch Two Sample t-test
data: log10(hum2$growth) and log10(lig2$growth)
t = 4.8606, df = 97.58, p-value = 4.473e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.4247076 1.0108416
sample estimates:
mean of x mean of y
7.000225 6.282451
t.test(log10(hum3$growth),log10(lig3$growth))
Welch Two Sample t-test
data: log10(hum3$growth) and log10(lig3$growth)
t = 3.4215, df = 83.534, p-value = 0.0009662
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.2080219 0.7855227
sample estimates:
mean of x mean of y
7.170545 6.673773
t.test(log10(hum4$growth),log10(lig4$growth))
Welch Two Sample t-test
data: log10(hum4$growth) and log10(lig4$growth)
t = 2.9997, df = 74.706, p-value = 0.00367
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1331595 0.6597929
sample estimates:
mean of x mean of y
7.343272 6.946796
t.test(log10(hum5$growth),log10(lig5$growth))
Welch Two Sample t-test
data: log10(hum5$growth) and log10(lig5$growth)
t = 3.1455, df = 56.335, p-value = 0.002647
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1801401 0.8117505
sample estimates:
mean of x mean of y
7.477756 6.981811
hum1 = out2[out2$fire_day ==1 & out2$eco1==11,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==11,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==11,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==11,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==11,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==11,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==11,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==11,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==11,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==11,]
mean(hum1$growth)
[1] 21604386
mean(hum2$growth)
[1] 43011405
mean(hum3$growth)
[1] 71836206
mean(hum4$growth)
[1] 110978456
mean(hum5$growth)
[1] 156565228
mean(lig1$growth, na.rm=T)
[1] 12710105
mean(lig2$growth, na.rm=T)
[1] 24454241
mean(lig3$growth, na.rm=T)
[1] 39894008
mean(lig4$growth, na.rm=T)
[1] NaN
mean(lig5$growth, na.rm=T)
[1] NaN
#t.test(log10(hum1$growth),log10(lig1$growth))
#t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
Error in t.test.default(log10(hum3$growth), log10(lig3$growth)) :
not enough 'y' observations
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 9)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 9))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 9)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 9))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 5 & month<10)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 5 & month<10))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
are ROS the same for light & human under the same conditions
summary(lm(out1$vpd~log(out1$ros_km+1),na.omit=T))
In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
extra argument ‘na.omit’ will be disregarded
Call:
lm(formula = out1$vpd ~ log(out1$ros_km + 1), na.omit = T)
Residuals:
Min 1Q Median 3Q Max
-2.11204 -0.40156 -0.05854 0.30543 2.12623
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.44554 0.02221 65.09 <2e-16 ***
log(out1$ros_km + 1) 0.52491 0.03383 15.52 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5888 on 1534 degrees of freedom
(265 observations deleted due to missingness)
Multiple R-squared: 0.1356, Adjusted R-squared: 0.1351
F-statistic: 240.7 on 1 and 1534 DF, p-value: < 2.2e-16
analysis of the first day
load data
daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/all_ignitions_V3.txt",header=T)
res=as.data.frame(daily_res)
res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$total_area =as.numeric(as.character(res$total_area))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))
res$biomass =as.numeric(as.character(res$biomass))
res = res[-1,]
res$human[res$cause ==1] =1
res$human[res$cause !=1 & res$cause !=14] =0
analysis
out1 = res[res$cause !=1 & res$cause != 14,]
out2 = res[res$cause ==1,]
out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1) #1=lightning; 14=unknown; 7=arson
out1 = subset(res,eco1 == 11& res$cause !=1 & res$cause != 14)
out2 = subset(res,eco1 == 11 & res$cause ==1)
out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14 & res$mont > 5 & res$mont < 10 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1& res$mont > 5 & res$mont < 10) #1=lightning; 14=unknown; 7=arson
t.test(out1$bi,out2$bi)
t.test(out1$erc,out2$erc)
t.test(out1$etr,out2$etr)
t.test(out1$fm100,out2$fm100)
t.test(out1$fm1000,out2$fm1000)
t.test(out1$pet,out2$pet)
t.test(out1$pr,out2$pr)
t.test(out1$rmax,out2$rmax)
t.test(out1$rmin,out2$rmin)
t.test(out1$th,out2$th)
t.test(out1$tmmn,out2$tmmn)
t.test(out1$tmmx,out2$tmmx)
t.test(out1$vpd,out2$vpd)
t.test(out1$vs,out2$vs)
t.test(out1$ws,out2$ws)
t.test(out1$biomass,out2$biomass)
t.test(out1$mean_land,out2$mean_land)
t.test(log10(out1$total_area),log10(out2$total_area))
ta = table(res$human,res$mont)
ps = barplot(ta, beside=TRUE, ylab="number of fires",xpd=T,xlab= "month", xaxt='n',ylim=c(0,300), axis.lty=1,cex.lab = 1.4,cex.axis = 1.2 )
axis(1,at=c(2,5,8,11,14,17,20,23,26,29,32,35), labels =c(1:12),xlim=c(0,36),xpd=F ,cex.lab = 1.4,cex.axis = 1.3)
legend("topright",c("human","lightning"),fill = c("grey", "black"), bty="n",cex=1.4)
text(1,285,"a)",cex=1.8)
out1 = subset(res,eco1 == 6 |eco1 == 7) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,eco1 == 11)
ta = table(out1$human,out1$mont)
ps = barplot(ta, beside=TRUE, ylab="number of fires",xpd=T, xaxt='n',xlab= "month", ylim=c(0,200), axis.lty=1,cex.lab = 1.4,cex.axis = 1.2 )
axis(1,at=c(2,5,8,11,14,17,20,23,26,29,32,35), labels =c(1:12),xlim=c(0,36),xpd=F,cex.lab = 1.4,cex.axis = 1.3 )
legend("topright",c("human","lightning"),fill = c("grey", "black"), bty="n",cex=1.4)
text(1,190,"b)",cex=1.8)
ta = table(out2$human,out2$mont)
ps = barplot(ta, beside=TRUE, ylab="number of fires",xpd=T,xaxt='n',xlab= "month", ylim=c(0,200), axis.lty=1,cex.lab = 1.4,cex.axis = 1.2 )
axis(1,at=c(2,5,8,11,14,17,20,23,26,29,32,35), labels =c(1:12),xlim=c(0,36),xpd=F,cex.lab = 1.4,cex.axis = 1.3 )
legend("topright",c("human","lightning"),fill = c("grey", "black"), bty="n",cex=1.4)
text(1,190,"c)",cex=1.8)
plot(data_forest$ros_km, data_forest$mean_BA_red,log="x",xlab="Rate-of-Spread (km/day)",ylab="Tree mortality (%)",xlim=c(0.005,30),xaxt="n",cex.axis=1.4 ,cex.lab=1.4,cex=0.8, col="black")
46 x values <= 0 omitted from logarithmic plot
axis(1,at=marks,labels=marks,cex.axis=1.4 )
points(data_test2$ros_km, data_test2$mean_BA_red,cex=0.8, col="darkgrey")
points(data_test1$ros_km, data_test1$mean_BA_red,cex=0.8, col="orange")
lines(lowess(data_forest$ros_km, data_forest$mean_BA_red, f=0.45),col="black", lwd=3)
lines(lowess(data_test1$ros_km, data_test1$mean_BA_red, f=0.45),col="darkgoldenrod3", lwd=3)
lines(lowess(data_test2$ros_km, data_test2$mean_BA_red, f=0.45),col="gray40", lwd=3)
NA
NA
quan = quantile(data_s1$ros_km,0.5)
fast = data_s1[data_s1$ros_km > quan,]
slow = data_s1[data_s1$ros_km < quan,]
fast_hum = fast[fast$human == 1,]
sum(fast$size)/sum(data_s1$size)
[1] 0.9084316
length(fast$size)/length(data_s1$size)
[1] 0.5
sum((data_s1$mean_BA_red*data_s1$size))/(sum(data_s1$size))
[1] 49.03327
mean(data_s1$mean_BA_red)
[1] 25.64874
sum(fast_hum$size, na.rm=T)/sum(fast$size)
[1] 0.4463899
length(fast_hum$size)/length(fast$size)
[1] 0.3658537
sum((fast$mean_BA_red*fast$size))/(sum(fast$size))
[1] 50.26964
mean(fast$mean_BA_red)
[1] 35.06328
sum((slow$mean_BA_red*slow$size))/(sum(slow$size))
[1] 36.76754
mean(slow$mean_BA_red)
[1] 16.2342
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
```